#include<stdio.h>
#include<math.h>
#include<gmp.h>
#include<time.h>

int main(void){
	mpf_t a0, b0, t0, p0;
	mpf_t an, bn, tn, pn, pi;	//double an, bn, tn, pn, pi;
	mpf_t auxA1, auxB1, auxT1, auxP1;	//double an, bn, tn, pn, pi;
	mpf_t auxA2, auxB2, auxT2, auxP2;	//double an, bn, tn, pn, pi;
	mpf_t pi_correct, precisao;
	
	int i = 0;
	FILE *fp;
	time_t t1, t2;
		
	fp = fopen("pi_input.txt", "r");
	if(fp == NULL){
		printf("Erro na abertura do arquivo!\n");
		return -1;
	}
	(void) time(&t1);
	/*inicializando as variáveis*/
	mpf_init(a0); 
	mpf_init(b0);
	mpf_init(t0); 
	mpf_init(p0);
	mpf_init(an); 
	mpf_init(bn);
	mpf_init(tn); 
	mpf_init(pn);
	mpf_init(pi);
	mpf_init(pi_correct);
	/*inicializando as variáveis auxilixares*/
	mpf_init(auxA1); 
	mpf_init(auxB1);
	mpf_init(auxT1); 
	mpf_init(auxP1);
	mpf_init(auxA2); 
	mpf_init(auxB2);
	mpf_init(auxT2); 
	mpf_init(auxP2);
	
	/*define a precisao*/
//	mpf_set_default_prec(33219280);
mpf_init_set_ui(precisao, 10);
mpf_pow_ui(precisao, precisao, 10000000); 
mpf_ui_div(precisao, 1, precisao);
	
	/*inicializando valores*/
	mpf_set_d(a0, 1.0); 
	mpf_set_d(b0, 0.707106781); 	
	mpf_set_d(t0, 0.25); 	
	mpf_set_d(p0, 1.0);
	
	printf("Retorno do input file = %d\n", mpf_inp_str (pi_correct, fp, 10)); 

	/*Impressão na tela*/
/*	printf("Printing initialized variables:\n");
	mpf_out_str (stdout, 10, 10, a0); printf ("\n");
	mpf_out_str (stdout, 10, 10, b0); printf ("\n");
	mpf_out_str (stdout, 10, 10, t0); printf ("\n");
	mpf_out_str (stdout, 10, 10, p0); printf ("\n\n");
	mpf_out_str (stdout, 10, 10, pi_correct); printf ("\n\n");
	*/
	
	do{
		/*an = (a0+b0)/2;*/
		mpf_add(auxA1, a0, b0);
		mpf_div_ui(auxA2, auxA1, 2.0);

		mpf_set(an, auxA2);

		/*bn = sqrt(a0*b0);*/
		mpf_mul (auxB1, a0, b0);
		mpf_sqrt(auxB2, auxB1);

		mpf_set(bn, auxB2);	

		/*tn = t0 - p0*(a0 - an)*(a0 - an);*/
		mpf_sub (auxT1, a0, an);	
		mpf_pow_ui (auxT2, auxT1, 2);
		mpf_mul (auxT1, p0, auxT2);
		mpf_sub (auxT2, t0, auxT1);	

		mpf_set(tn, auxT2);	

		/*pn = 2 * p0;*/
		mpf_mul_ui (pn, p0, 2);
		
		/*pi = (an + bn)*(an + bn)/(4*tn);*/	
		mpf_add(auxA1, an, bn);
		mpf_pow_ui (auxP1, auxA1, 2);
		mpf_mul_ui (auxP2, tn, 4);
		mpf_div (pi, auxP1, auxP2);

		mpf_set(a0, an);
		mpf_set(b0, bn);
		mpf_set(t0, tn);
		mpf_set(p0, pn);
		
		//printf("Pi [iteracao == %d] = ", i);
		//mpf_out_str (stdout, 10, 10, pi); printf ("\n\n");
		
		mpf_sub(auxA1, pi, pi_correct); 
		mpf_abs(auxA1, auxA1);
		i++;
	} while(mpf_cmp (auxA1, precisao) <= 0 || i < 159);

	(void) time(&t2);
	
	printf("Resultado final: Pi = "); mpf_out_str (stdout, 10, 10, pi); printf ("\n");
	printf ("Numero de iteracoes = %d\nTempo de calculo = %f segundos\n", i, difftime(t2,t1));

	return 0;
}
